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Abstract 



This paper studies the effects of distributed delay coupling on the dynamics in a system 
of non-identical coupled Stuart-Landau oscillators. For uniform and gamma delay distribution 
^ ' kernels, conditions for amplitude death are obtained in terms of average frequency, frequency 

detuning and parameters of the coupling, including coupling strength and phase, as well as the 
mean time delay and the width of the delay distribution. To gain further insight into the dynam- 
ics inside amplitude death regions, eigenvalues of the corresponding characteristic equations are 
computed numerically. Oscillatory dynamics of the system is also investigated using amplitude 
and phase representation. Various branches of phase-locked solutions are identified, and their 
stability is analysed for different types of delay distributions. 



1 Introduction 



Dynamics of many complex physical, biological, and engineering systems can be effectively modelled 
mathematically using ensembles of coupled oscillators (Kuramoto 1984, Pikovsky et al. 2001). A 
significant advantage of such an approach over other methodologies lies in the possibility to identify 
and study regimes with particular types of behaviour, such as emergence and stability of cooperation 
and (de)synchronization, different kinds of spatial patterns, travelling waves etc. Such analyses 
often provide very important insights for practical applications. Our understanding and treatment 
of various brain pathologies and deficiencies, such as Parkinson's, Alzheimer's, epilepsy, heavily 
relies on the analysis of synchronization of oscillating neural populations (Uhlhaas & Singer 2006, 
Popovych et al. 2005, Schnitzler & Gross 2005). Recent work on laser communication networks 
has extensively used coupled oscillator models to study the dynamics of in-phase or anti-phase 
and complete chaos synchronization (Heil et al. 2001, Flunkert et al. 2009, Hicke et al. 2011, 
Flunkert & Scholl 2012, Kozyreff et al. 2000). Chimera states, where a network of oscillators 
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splits into coexisting domains of coherent, phase-locked and incoherent, desynchronized behaviour, 
have also been observed (Abrams & Strogatz 2004, Hagerstrom et al. 2012, Omelchenko et al. 
2011, Kuramoto & Battogtokh 2002, Tinsley et al. 2012). These studies have shown that coupling 
between different elements can play a dual role in the dynamics: it can both lead to suppression of 
oscillations, and it can also facilitate a certain degree of synchronization between different elements. 

When analysing the dynamics of coupled systems, it is important to take into account the 
fact that in many cases the coupling between different elements is not instantaneous. Time delays 
associated with such coupling can have a major impact on the overall dynamics of the system. In 
practical examples these time delays are often associated with delays in propagation or processing 
of signals, response time of mechanical actuators, translation and transcription time in genetic 
oscillators (Atay 2010, Just et al. 2010, Kyrychko & Hogan 2010, Scholl & Schuster 2008). In 
all these examples, explicit inclusion of time delays in the model has provided a more realistic 
and accurate representation of the system under consideration. From a mathematical perspective, 
differential equations with time delays are infinite-dimensional dynamical systems, which makes 
their analysis, as well as numerical simulations, quite involved (Diekmann et al. 1995, Erneux 
2009). 

In the case when the coupling between oscillators is sufficiently weak, it is possible to reduce the 
full system to a phase-only model (Kuramoto 1984). Significant amount of research has been done 
over the years on the analysis of phase oscillators with different kinds of time-delayed coupling 
(D'Huys et al. 2008, Earl & Strogatz 2003, Erneux & Glorieux 2010, Nakamura 1994, Niebur 
1991, Pikovsky et al. 2001, Popovych et al. 2010, Schuster & Wagner 1989, Sen et al. 2010). 
Such systems of phase oscillators demonstrate a rich variety of dynamical regimes, including chaos, 
synchronization, splay and chimera states, characterized by co-existence of coherent and incoherent 
states. In many cases, coupling in the phase models leads to phase entrainment and the emergence 
of so-called phase-locked solutions, where all oscillators start to oscillate with the same common 
frequency and have a constant shift between their phases. Existence and stability of such phase- 
locked solutions have been studied in a number of systems of oscillators with different kinds of local 
and non-local coupling (Kuramoto & Battogtokh 2002, Sen et al. 2010). 

For sufficiently strong coupling between oscillators, one can observe another interesting and 
important aspect of dynamics of coupled oscillator systems: the ability of external coupling to 
suppress otherwise stable periodic oscillations. This was first discovered by Bar-Eli in the context 
of chemical reactions with instantaneous coupling (Bar-Eli 1985), and it has been shown for a 
system of two coupled Stuart-Landau oscillators that the 'Bar-Eli effect' can only occur, provided 
the coupling strength and the frequency detuning of the two oscillators are both sufficiently large 
(Aronson et al. 1990, Ermentrout 1990, Mirollo &: Strogatz 1990). Later it was discovered that when 
the coupling is time-delayed, it is possible to achieve suppression of oscillations and stabilization 
of an unstable fixed point even for identical oscillators (Ramana Reddy et al. 1998,1999). This 
phenomenon, named amplitude death (Ramana Reddy et al. 1998), oscillator death, or 'death 
by delay' (Strogatz 1998), has been demonstrated experimentally in nonlinear electronic circuits 
(Ramana Reddy 2000), in the dynamics of slime mold Physarum polycephalum (Herrero et al. 
2000) and thermo-optical oscillators linearly coupled by heat transfer (Takamatsu et al. 2000). 
Amplitude death has been subsequently studied for a number of different systems and couplings 
(Ramana Dodla et al. 2004, Karnatak et al. 2010, Sen et al. 2005, 2010, Choe et al. 2007, Scholl 
et al. 2009, Zuo et al. 2011, 2012, Bar et al. 2012). 

Despite successes in the analysis of systems of coupled oscillators with time-delayed coupling, one 
of the limitations of the majority of this research has been the restriction on the type of time-delayed 
coupling, which is usually taken to be in the form of one or several constant time delays. At the same 
time, in many realistic systems time delays themselves are not constant (Gjurchinovski &: Urumov 
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2008, 2010) and may either vary depending on the values of system variables (state-dependent 
delays) or just not be explicitly known. In order to account for such situations mathematically, one 
can use the formalism of distributed time delays, where the time delay is represented through an 
integral kernel describing a particular delay distribution (Bernard 2001, Campbell & Jessop 2009). 
Distributed time delay has been successfully used to describe situations when only an approximate 
value of time delay is known in engineering experiments (Kiss &; Krauskopf 2011, Michiels et al. 
2005), for modelling distributions of waiting times in epidemiological models (Blyuss & Kyrychko 
2010), maturation periods in population and ecological models (Faria 2010, Gourley 2003), as well 
as in model of traffic dynamics (Sipahi et al. 2008), neural systems (Eurich et al. 2005), predator- 
prey and food webs (Thiel et al. 2003). 

In this paper we investigate amplitude death and phase dynamics in a system of coupled Stuart- 
Landau oscillators with distributed delay coupling. This system is a prototype for dynamics near 
a supercritical Hopf bifurcation, and in this capacity it captures essential features of many realistic 
systems in such a regime. The corresponding mathematical model can be written in the form 



Zl {t) = (1 + iwi)zi(i) - \zi{t)\ 2 Zl {t) + Ke ie 



o 



Z2(t) = (1 + ico 2 )z 2 (t) - \z 2 {t)\ 2 z 2 {t) + Ke id 



g{t')z 2 {t-t')dt' - Zl {t) 



g{t')z x {t - t')dt' - z 2 (t) 



(1) 



o 



where Z\ 2 G C, w\ 2 are frequencies of two oscillators, K E M + and € M are the strength and the 
phase of coupling, respectively, and g(-) is a kernel of delay distribution. The kernel g is taken to 
be positive-definite and normalized to unity: 

/•oo 

g(u) > 0, / g(u)du = 1. 
J o 

The case g(u) = S(u) corresponds to instantaneous coupling (z 2 — z\) describing a situation when 
oscillators interact without any delay. For g(u) = 5(u — t), the coupling takes the form of a discrete 
time delay: [z 2 (t — r) — z%(t)]. Atay (2003) has analysed this system for the case of zero coupling 
phase and a uniform delay distribution kernel analytically for the case of identical oscillators and 
numerically for non-identical oscillators, and identified regimes of amplitude death in terms of 
coupling strength and mean time delay. Kyrychko et al. (2011) have analysed amplitude death in 
model ([1]) with lv\ = oj 2 = ujq for uniform and gamma distributed delay kernels and non-zero phase. 

The outline of this paper is as follows. In the next section we analyse stability of the triv- 
ial equlibrium of system (JTJ) and find regions of amplitude death depending on parameters of the 
coupling and delay distribution. The eigenvalues of the corresponding characteristic equations de- 
termining the stability of the steady state are computed numerically to gain a better understanding 
of system dynamics inside amplitude death regimes. Section 3 contains analysis of phase-locked 
solutions of system ([T|) for uniform and gamma distributions in the case of constant equal ampli- 
tudes of oscillations (when the system reduces to a coupled system of Kuramoto oscillators with 
distributed-delay coupling), as well as in a general case of arbitrary amplitudes. In each case, we 
identify branches of phase-locked solutions and numerically compute their stability. The paper 
concludes with a discussion of results as well as possible further developments of this work. 
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2 Amplitude death 



To study the possibility of amplitude death in the system ([T]), we linearize this system near the 
trivial steady state Z12 = 0- The corresponding characteristic equation is given by 

1 +iui - Ke ie " A) (l + iu 2 - Ke ie - a) - K 2 e 2W {{Cg}(\)\ 2 = 0, (2) 
where A is an eigenvalue of the Jacobian, and 

{Cg}(s) = / e- su g(u)du, (3) 



is the Laplace transform of the function g{u). Atay (2003) has investigated the case 8 = an- 
alytically for the case of identical oscillators oj\ = uj 2 = ^o> an d numerically for oj\ ^ io 2 and a 
uniform delay distribution kernel. More recently, Kyrychko et al. (2011) have studied analytically 
and numerically the case of ui\ = io 2 = ujo an d a non- vanishing coupling phase 9^0 with uniform 
and gamma delay distributions. 

To make further analytical progress, it is instructive to specify a particular choice of the delay 
kernel. As a first example, we consider a uniformly distributed kernel 

— for r — p < u < t + p, 
.</(») { 9 ' (4) 

elsewhere. 

This distribution has the mean time delay 



/>oo 

=< r >= / ug{u)du = r, 
Jo 



and the variance 

f°° P 2 
a 2 = I (u- T m ) 2 g(u)du = —. (5) 

Jo 

In the case of uniformly distributed kernel (j3J), it is quite easy to compute the Laplace transform 
of the distribution g(u) as: 

1 ^-Ar (\p _ c -\p\ _ ^-Ar sinh ( A / 5 ) 



and this transforms the characteristic equation (|2|) into 

(l + kj x - Ke w - A) (l + iuj 2 - Ke w - a) = K 2 e 2ld e - 2XT 



sinh(Ap)" 1 *" 
Ap . 



(6) 



Since the roots of the characteristic equation (|6|) are complex- valued, stability of the trivial steady 
state can only change if some of these eigenvalues cross the imaginary axis. To this end, we can 
look for characteristic roots in the form A = iu. Substituting this into the characteristic equation 
([6]) and separating real and imaginary parts gives the following system of equations for (K, r) as 
parameterized by the Hopf frequency w. 

A 2 

(1 -KcosO) 2 - (57 - u - KsmO) 2 + — = K 2 5( P ,uj) cos[2(fl - urr)], 

4 (7) 

2(1 - K cos«) (57 - K sin6» - w) = K 2 5(p, u) sin[2(6» - ut)], 
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Figure 1: (Colour online) Regimes of amplitude death depending on the coupling strength K and 
r for the uniform distribution kernel with 9 = 0, p = (discrete time delay), and U = 10. Colour 
code denotes [-max{Re(A)}] for max{Re(A)} < 0. (a) A = 0. (b) A = 1.45. (c) A = 1.65. (d) 
A = 1.75. 



where 



S(p,u>) 



sin(up) 
up 



and we have introduced the frequency mismatch (detuning) A and the mean frequency oo as 

_ UJi + UJ 2 

A = wi-o;2, w = 5 • 

Squaring and adding the two equations in ([7]) gives a single quartic equation for the coupling 
strength K 

A^l 2 

4 J (8) 



(1 - K cos 9f - (w - oj - K sin Of + 



+4(1 - K cos 6) 2 {uo -Ksm6- w) 2 = K 4 5 2 ( P ,uj). 

In a similar way, dividing the two equations in (J7|) gives the equation for the time delay r at the 
Hopf bifurcation as 

2(1 - Kcos6)(u- oj - KsinO) 



tan[2(0 - ur)] 



(1 - KcosO) 2 - (5J - w- ifsine) 2 + A 2 /4' 



(9) 
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Figure 2: (Colour online) Amplitude death region = 0, r = 0.15, lo = 10, is bounded by the 
two surfaces in figure (a). Lines in (b) show cross-sections for A = 3 (solid), A = 5 (dashed), and 
A = 7 (dash-dotted). 



When the coupling phase vanishes (0 = 0), the above equations for (K,r) simplify to 



A 2 

1 -K) 2 - (lo-lo) 2 + — 



+ 4(1 - K) 2 (W - uf = K*S\p, lo), 



4r2, 



tan(2wr) 



2(1 - K){lo-lo) 



(1-K) 2 - (ld-Lo) 2 + A 2 /4' 



For A = 0, we have two identical oscillators with the same frequency lo\ = L02 = 
this into equation (|8|) gives the equation for the Hopf frequency in the form 

(1 - K cos Of + (5J - w - K sin 9) 2 = K 2 5{p, u) , 

and the equation for the critical time delay r at the Hopf bifurcation reduces to 

2(l-Kcos9)(oJ-u)-Ksm6) 2z 
(1-K cos 8) 2 - (57 - u - K sin 9) 2 ~ 1 - z 2 ' 



(10) 

ojq. Substituting 
(11) 



tan[2(0 - lot)] 



u — oj — K sin 9 
Z ~ I- K cos 6 

Using the trigonometric identity tan 2a = 2tana/(l — tan 2 a), we find 

cJ — to — K sm8 

tan(6/ — lot) - 



1 - K cos 9 



(12) 



Equations (llip and (I12p have been recently studied by Kyrychko et al. (2011), where the effects of 
the width of delay distribution p, as well as coupling strength K and the coupling phase 9 on the 
amplitude death were investigated. In the case of vanishing coupling phase (0 = 0), these equations 
reduce even further to a system studied by Atay (2003). 
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Figure 3: (Colour online) Regimes of amplitude death depending on the coupling strength K and 
coupling phase for the uniform delay distribution kernel with r = 0.15, p = 0.002 and U = 10. 
Colour code denotes [-max{Re(A)}] for max{Re(A)} < 0. (a) A = 0. (b) A = 3. (c) A = 6. (d) 
A = 9. 

To illustrate the effects of varying the coupling strength K and the time delay r on the (in)stability 
of the trivial steady state, we now compute the stability boundaries dHJ)-© as parameterized by the 
Hopf frequency u. Besides the stability boundaries themselves, which enclose the amplitude death 
regions, we also compute the maximum real part of the eigenvalues using the traceDDE package 
in Matlab. In order to compute these eigenvalues, we introduce real variables z\ r i and z% r i, where 
z\ = Zi r + izu and z<i = zi r + i&n, and rewrite the linearized system (SL) with the distributed 
kernel Q as 

z(t) = L z(t) + — / Mz(t + s)ds, (13) 

lp J-(r+p) 



where 



Ln 



Ni 2 
2 N 2 
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Figure 4: (Colour online) Regimes of amplitude death depending on the coupling strength K and 
a for the weak delay distribution kernel (fT5|) with = 0, p = 1 and tJ = 10. Colour code denotes 
[-max{Re(A)}] for max{Re(A)} < 0. (a) A = 0. (b) A = 1.45. (c) A = 1.65. (d) A = 1.75. 



and also 



1 - 



M : 

if cos (9 
- K sin 6 



2 R 
R 2 



K sin 6 — uji 
1 - KcosO 



R 



cos y 
— sin I 



sin t 

COS I 



No 



1 - Kcos9 
UJ2 — K sin 6 



if sin 9 — L02 
1 - if cos 9 



and 2 denotes a 2 x 2 zero matrix. When p = 0, the last term in the system (|13j) turns into 
KMz(t — r), which describes the system with a single discrete time delay r. System (|13p is in the 
form in which it is amenable to the algorithms described in Breda et al. (2006) and implemented 
in traceDDE. 

First of all, we consider the case p = and 8 = 0, corresponding to a discrete time delay and 
vanishing coupling phase, as shown in Fig. [H This figure shows that as the frequency detuning 
increases, this leads to emergence of new islands of amplitude death, and for sufficiently high 
detuning A, these islands merge into a single continuous region in the plane of coupling strength K 
and average time delay r, where amplitude death can occur for an arbitrary value of r, provided K 
lies in the appropriate range. One can note that unlike the case of identical oscillators considered 
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Figure 5: (Colour online) Regimes of amplitude death depending on the coupling strength K and 
the coupling phase 9 for the weak delay distribution kernel (|15p with a = 10, p = 1 and uJ = 10. 
Colour code denotes [-max{Re(A)}] for max{Re(A)} < 0. (a) A = 0. (b) A = 2. (c) A = 5. (d) 
A = 9. 



in Kyrychko et al. (2011) in this situation the values of K needed to achieve stabilization of the 
trivial steady state for any r are in the lower part of the overall range of K values. Increase in the 
width of uniform delay distribution p leads to a corresponding increase in the region of amplitude 
death, as illustrated in Fig. [2j When we consider the impact of coupling phase on amplitude death, 
as shown in Fig. El it becomes clear that the largest range of admissible K values is attained for 
9 = 0, but overall the range of coupling phases, for which amplitude death is possible, increases 
with the frequency detuning A. It is worth noting, however, that unlike the situation considered 
in Kyrychko et al. (2011) where increase in the width of distribution p led to the asymmetry 
of amplitude death regions with regards to the coupling phase, as A is increased, the region of 
amplitude death remains symmetric in 9 (for small p). 

The second example we consider is that of a gamma distribution, which can be written as 



9(u) 



T{p) 



(14) 



with a,p > 0, and T(p) being the Euler gamma function defined by T(0) = 1 and T(p + 1) = pT{p). 
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Figure 6: (Colour online) Regimes of amplitude death depending on the coupling strength K and a 
for the strong delay distribution kernel (fT5j) with 6 = 0, p = 2 and U = 10. (a) Stability boundary, 
(b) Stability boundary for A = 2. Colour code denotes [— max{Re(A)}] for max{Re(A)} < 0. 



For integer powers p, this can be equivalently written as 

9{u) 



(p-1)! 



(15) 



For p = 1 this is simply an exponential distribution (also called a weak delay kernel) with the 
maximum contribution to the coupling coming from the present values of variables z\ and z<i- For 
p > 1 (known as strong delay kernel in the case p = 2), the biggest influence on the coupling at any 
moment of time t is from the values of at t — (p — l)/a. The delay distribution (|15p has the 
mean time delay 

p 



ug(u)du = — , 
o « 



(16) 



and the variance 



a 



(u - T m )' 2 g(u)du = —. 



a- 



When studying stability of the trivial steady state of the system ([T]) with the delay distribution 
kernel (fT5j) . one could use the same strategy as the one described for uniform distribution. The 
Laplace transform of the distribution kernel in this case is given by 



{A?}(A) 



a 1 



(\ + a)P' 



Substituting this into the characteristic equation ([2]) yields a polynomial equation for A: 



1 + iuj\ — Ke ie — A ) ( 1 + iw 2 - Ke M - A ) (A + a) Zp = K l e llt> a 



\2p 



-2 2i6 2p 



(17) 



In the Appendix we show how the same characteristic equation can be derived using a linear chain 
trick without resorting to the Laplace transform. 

Figure d] illustrates the regions of amplitude death for a weak delay distribution kernel (|15p 
with p = 1, vanishing coupling phase and increasing frequency detuning A. Similarly to the case 



10 



10 



K 



o 




10 



K 



o 







10 r 










0.8 


K 

5 - 




0.4 













-71/2 -71/4 TU/4 71/2 

e 




-n/2 -7E/4 m 7E/2 

e 




-7U/2 -7U/4 71/4 7C/2 

6 







10 r 




1.2 






0.8 


K 

5 - 




0.4 














-7U/2 -71/4 7E/4 7C/2 

e 



Figure 7: (Colour online) Regimes of amplitude death depending on the coupling strength K and 
the coupling phase 9 for the strong delay distribution kernel (I15p with a = 1, p = 2 and W = 10. (a) 
A = 0. (b) A = 2. (c) A = 4. (d) A = 7. Colour code denotes [-max{Re(A)}] for max{Re(A)} < 0. 



of uniform delay distribution, as A increases, new islands of amplitude death appear, and eventually 
they merge into a single continuos region of amplitude death. Since in this figure, the regions of 
amplitude death are plotted in terms of a, which according to (fT6j) is the inverse average time 
delay r, this implies that for the case of a single connected region of amplitude death in the (a, K) 
plane, amplitude death can happen for an arbitrarily small value of the average time delay, provided 
frequency detuning is sufficiently large. In Fig. [5] we illustrate how amplitude death regions depend 
on the coupling phase. This figure suggests that for the same average time delay, provided the 
coupling phase is sufficiently negative, it is possible to achieve amplitude death for an arbitrary 
value of the coupling strength K starting from some minimal value. The regions of amplitude death 
are strongly asymmetric in 9, and amplitude death is possible for higher positive values of 9 for 
larger frequency detuning A. 

When one considers a strong distribution kernel (|15p with p = 2 with vanishing coupling phase, 
there is a minimum value of the average time delay required to achieve amplitude death, as shown 
in Fig. [6j This figure also suggests that the actual value of the minimum average time delay 
increases with the increasing frequency detuning A. Figure [7] illustrates how the coupling phase 
affects amplitude death. One can notice that similar to the case of a weak distribution kernel, the 
regions of amplitude death are asymmetric in 8. However, there are two major differences from 
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the case p = 1, namely, that the regions of amplitude death are closed in the (K,9) parameter 
plane, and these regions shrink with increasing A rather than grow as was the case for p = 1. The 
implication is that now it is not possible to find a coupling phase, for which amplitude death could 
be achieved for an arbitrary value of the coupling strength K. 



3 Phase-locked solutions 



In the previous section we analysed the situation when the distributed time delay in the coupling 
leads to the destruction of stable limit cycles and stabilization of the previously unstable fixed 
point. In order to understand the phase dynamics of the system, we introduce new real variables 
ORi(i),ft 2 (t),<Mi),<M*)) as follows 

*i(t) = Ri{t)e^\ z 2 (t) = R 2 {t)e i ^ t \ 

Here R\2 > and (f>i t 2 denote the amplitudes and phases of the two oscillators, respectively. 
Substituting this representation into the system ([I]) yields the following system of equations for the 
amplitude and phase variables 

r poo 

R 1 = (l-Rl)R 1 + K / g(t')R 2 (t - t') cos[(t> 2 (t - t') - ((>i(t) + 6]dt' -Rx cos 6 , 

Jo 

r r°° 

R 2 = (l - Rj) R 2 + K \ g(t')R!(t - t') cos[(t>i(t - t') - (f) 2 (t) + 0]dt' -R 2 cos 

Jo 

r r°° 

Rij>i = RiUi+K / g(t')R 2 (t-t')sm[^ 2 (t-t') - <f>i(t) + 0}dt' - R x sin 



(18) 



R 2 <t>2 = R 2 u 2 + K 



g{t')R 1 {t-t')sm[(t) 1 {t-t') - (j) 2 {t) + e\dt' - R 2 sin 6 



Significant amount of work has been done on the analysis of this system for the case of instantaneous 
coupling g(s) = 5(s) or discrete time delay g(s) = 5(s —r). In both of these cases, it has been 
shown that besides amplitude death, the system can exhibit a number of other interesting solutions, 
such as phase-locked (also known as frequency-locked) and phase drift solutions. In this section we 
consider effects of distributed-delay coupling on the dynamics of such solutions, which have so far 
remained unexplored. 



3.1 Constant amplitude dynamics 

As a first step in the analysis of system f)18[) . we assume that the amplitudes of both oscillators are 
equal to each other and constant: R\{t) = R 2 (t) = const for all times. In this case, neglecting am- 
plitude dynamics, the system (fT8|) reduces to a system of two Kuramoto oscillators with distributed 
delay coupling 



uix + K 



0J2 + K 



g(t')sm[(j)2(t - t') - 4>x(t) + 9}dt' - sin I 
3(0 sin [0i (i _ t') - 4> 2 (t) + 9]dt' - sin I 



(19) 
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It is instructive to consider this system using the variables of the phase difference tp = 4>i — <p2 and 
the mean phase tp = {4>\ + $2)/ 2: 



(p = uj + K 



' g(t') sin (<p(t - t') - <p(t) + 9) cos t^l±^LJl dt > _ sin , 



i> = A - 2K I g(t') cos (<p(t - - <p(t) + 9) sin *£M1±M l !l dt > . 



(20) 



The last equation suggests that when A > 2K, this system exhibits phase drift (i.e. unbounded 
growth of the difference between the phases of two oscillators) independently of the delay distribu- 
tion kernel, and hence phase locking can only occur for A < 2K. Following Schuster and Wagner 
(1989), we look for solutions of the system (|20p in the form 



{tp\il)*) = (nt + const, /3), 



(21) 



where O is the ensemble frequency of oscillations, and /3 is the constant phase shift between the 
two oscillators. Substituting this into the system ([20]) yields the value of the phase shift as 



arcsm 



A 

2if 



7r — arcsm 



2K 



F c (-n,9)- 1 



(22) 



where we have introduced the notation 

poo poo 

F c (a,b) = / g(t') cos(at' + b)dt' , F s (a,b) = / g(t') s'm(at' + b)dt' . 
Jo Jo 

The new common frequency Q satisfies the transcendental equation 

F s (-n,0) 



f±(Q) =57- ft- Ksm0± 



„K 2 F?(-n,6)- — 

F c (-n,0) V cV ' 4 



0, 



(23) 



where the plus and minus sign correspond to the first and second value of j3 in (|22|) . respectively. 
It follows from this equation that possible solutions for ft can only lie in the admissible range 



F c (-n,0f > 



A 2 
4K 2 ' 

Stability of the phase-locked solution (|21[) is determined by the roots of the corresponding charac- 
teristic equation 



A 2 + 2K cos PF C (-Sl, 0) + K 2 [F c (-n, 9 - /3)F c (-ft, 9 + (3) 

-F c L (-ft, - p, A)F c (-ft, 9 + 0, A)] = 0, 



(24) 



where 



poo 

F^(a,b,z)= g{t') cos(at' + b) 
Jo 



e- zt 'dt', 
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Figure 8: (a) Function /±(^) from (j23[) for uniform distribution kernel (|4]) with a; = 10, r = 4, 
p = 0.1, = 0, A = 1, and K = 0.8 (solid), K = 1.2 (dashed), K = 1.5 (dash-dotted), (b) Function 
f±(£l) from (|23p for uniform distribution kernel (JH) with ZJ = 10, t = 4, p = 0.2, 9 = 0, K = 1.5 and 
A = (solid), A = 0.4 (dashed), A = 0.8 (dash-dotted), (c-f) Branches of phase-locked solutions 
(j2T|) for uniform delay kernel with 57 = 10, p = 0.2, 6 = 0, r = 4, and A = (c), A = 0.4 (d), 
A = 0.8 (e), A = 1.2. Solid lines denote stable branches, dashed lines denote unstable branches. 

and the superscript L refers to this integral representing the Laplace transform of a modified kernel 
g(s) cos(as + b). 

To make further progress, we have to use specific form of a distributed kernel. For the uniform 
distribution @, functions F c and F s can be computed as follows 

F c (a,b) = — sin (ap) cos (ar + b), F s (a,b) = — sin( a p) sin( ar + b), 
pa pa 

and the case of a discrete time delay can be recovered by setting p = 0. For the weak distribution 
kernel (|15p with p = 1, we have 

a cos b — a sin b . . a sin b + a cos 6 

F c (a, b) = a =- — , F s (a,b) = a- 



a 2 + a 2 



a 2 + a 2 
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Figure 9: (a) Function /±(f2) from f|23|) for the weak distribution kernel (|15p with p = 1, u = 10, 
a = l, = 0, A = 1 and K = 10 (solid), K = 20 (dashed), if = 50 (dash-dotted), (b) Function 
/±(fi) from ([23]) for the weak distribution kernel (JT5D with p = 1, 57 = 10, a = 1, 6> = 0, K = 200 
and A = (solid), A = 0.4 (dashed), A = 0.8 (dash-dotted), (c) Branches of phase-locked solutions 
(|2ip for the weak delay distribution kernel (|15p with tJ = 10, p = 1, and a = 1. All branches are 
unstable. 



and similarly, for the strong delay kernel (|15p with p = 2: 

2 (a 2 — a 2 ) cos 6 — 2crasin b 2 (° 2 ~~ ° 2 ) srn ^ + 2aa cos 6 

Fc(a ' 6) = a ^r^? ' F ^ = a wr^y • 

Figure [8] illustrates how the function f±(£l), whose roots determine the values of the collective 
frequency Q of phase-locked branches, depends on frequency detuning A and the coupling strength 
K . Computation of the stability of the branches shown in this Figure indicates that as the frequency 
detuning increases, branches of phase-locked solutions start to appear for higher values of the 
coupling strength. 
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Figure 10: (a) Function /±(0) from (|23|) for the strong distribution kernel (|15p with p = 2, u = 10, 
a = l, 6» = 0, A = 3 and K = 10 (solid), K = 20 (dashed), K = 50 (dash-dotted), (b) Function 
/±(^) from (|23p for the strong distribution kernel ()15p with p = 2, uJ = 10, a = 1, 6 = 0, .K' = 10 
and A = (solid), A = 4 (dashed), A = 8 (dash-dotted), (c-f) Branches of phase-locked solutions 
(f2Tj) for the strong delay distribution kernel (fT5j) with uJ = 10, p = 2, ct = 1, 6 = 0, and A = 
(c), A = 3 (d), A = 6 (e), A = 9. Solid lines denote stable branches, dashed lines denote unstable 
branches. 

Similar computation for the weak delay distribution kernel, shown in Fig. [91 suggests that in 
this case all branches of phase- locked solutions are unstable for any value of A. In the case of the 
strong distribution kernel, which is illustrated in Fig. [TPl the branches of phase-locked solutions 
are stable for a very narrow range of K near the lowest tip of the branch for sufficiently small 
detuning (see insets). As the detuning increases, this makes the branches of phase- locked solutions 
to appear for higher values of the coupling strength in a manner similar to the case of uniform 
delay distribution kernel, and all these branches are unstable. 

3.2 General phase-locked solutions 

In a more general case, when R\ and i?2 are not required to be equal to each other, one can still 
look for phase-locked solutions of the system (|18h in the form 

{R x {t),R 2 {t), 0i (t), <h{t)) = (Ri,R2, nt + p/2, nt - p/2), 

where R\2 are unknown constants, O is the new common frequency of oscillations, and /3 is the 
phase shift between the two oscillators. These solutions can be found as the roots of the following 
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Figure 11: Phase-locked solutions (J25J) for the uniform delay distribution kernel with p = 0.2, 
u = 10, r = 4, 6 = 0. (a) A = 0. (b) A = 0.2 (c), A = 0.4 (d), A = 0.6. Circles denote the branch 
with Ri = R 2 , lines denote the branches with i?i ^ R 2 . Solid lines (closed circles) denote stable 
branches, dashed lines (open circles) denote unstable branches. 



R 2 I g (f ') cos(- P -fit' + 6) dt' - R x cos 1 



Rx I g(t')cos(p-nt' + 6)dt'-R 2 cos6 



0, 



(25) 



system of equations 

Ri + K 

(1 - r£) r 2 + k 

RiQ = Riui + K 
R 2 tt = R 2 oj 2 + K 



It is easy to check that the constant amplitude phase-locked solutions f)21 1) are only the solutions 
of the full system ()25[) for the case of identical oscillators uj\ = uj 2 , which implies A = and the 
phase shift /3 of either zero or ir (describing in-phase or anti-phase solutions, respectively). 

Although it is not possible to find solutions of the system (|25p analytically, the phase-locked 
solutions can be computed numerically for each particular choice of delay kernel. Figure [TT] shows 



R 2 I g (tf ')sin(- P - Qtf + 0) dt' - R t sin 6 


Ri I g(t') sin(/3 - Sit' + 9) dt' - R 2 sin 6 
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Figure 12: Phase-locked solutions (|25[) for the weak delay distribution kernel (|15p with a = 1, 
p = 1, u = 10, e = 0. (a) A = 0. (b) A = 0.1 (c), A = 0.2 (d), A = 0.4. Circles denots the branch 
with R\ = i?2) black lines denote the branches with R\ ^ R 2 - Solid lines (closed circles) denote 
stable branches, dashed lines (open circles) denote unstable branches. 

the branches of phase-locked solutions (I25p for the uniform delay distribution kernel. For A = 0, 
we identify both the branch of equal amplitude phase-locked solution, and also an unstable branch 
with unequal amplitudes. As the frequency detuning increases, stable branches of phase-locked 
solutions are observed for higher values of the coupling strength. For all values of detuning, as K 
tends to zero, the ensemble frequency of branches of phase-locked solutions approach the values of 
Q = wi,2) which should be expected from the fact that for K = the system (I25p admits solutions 
(R!,R2,n) = (1,0, wi) and (R 1 ,R 2 ,U) = (0,1, cj 2 ). 

Figure H2] illustrates the branches of phase-locked solutions for the weak distribution kernel. 
Similarly to the case of uniform delay distribution kernel, for A = we find the branches of 
solutions with both equal and different amplitudes, and both of these branches are unstable. For 
higher values of frequency detuning, we have two unstable branches of phase-locked solutions located 
symmetrically around the average frequency uJ. The distance between these branches in terms of 
their ensemble frequency is increasing with increasing A. For the strong distribution kernel, the 
situation is similar, as shown in Fig. [I3j Once again, we observe two unstable branches of phase- 
locked solutions, whose frequencies are centred around uJ, and the distance between the branches 
grows with A. 
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Figure 13: Phase-locked solutions (|25p for the strong delay distribution kernel (|15p with a = 1, 
p = 2, oJ = 10, 6 = 0. (a) A = 0. (b) A = 0.1 (c), A = 0.2 (d), A = 0.4. Circles denote the branch 
with R\ = i?2) black lines denote the branches with R\ 7^ i?2- Solid lines (closed circles) denote 
stable branches, dashed lines (open circles) denote unstable branches. 



4 Discussion 

In this paper we have analysed amplitude death, as well as existence and stability of phase-locked 
oscillatory solutions in a generic model of Stuart-Landau oscillators with delay-distributed coupling. 
Regions of amplitude death have been identified for uniform and gamma distributions depending 
on the intrinsic system parameters (average frequency and detuning), as well as coupling strength, 
phase and the distribution parameters. Computation of boundaries of amplitude death regions, 
as well as stability eigenvalues for the trivial steady state, suggests that the more disparate the 
oscillators are (i.e., the higher the frequency detuning is), the larger are the regions of amplitude 
death in the case of uniform and weak delay distributions. However, for the strong delay distribution 
kernel, higher detuning corresponds to a smaller region of amplitude death. Special attention has 
been paid to the analysis of the effects of the coupling phase on amplitude death. While for uniform 
delay distribution, the maximum range of coupling strengths giving amplitude death is achieved 
for zero phase, in the case of gamma distribution, the non-zero coupling phase increases the range 
of such coupling strengths. 

To understand the dynamics beyond amplitude death, we have analysed the appearance and 
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stability of phase-locked solutions, characterized by both oscillators performing oscillations with 
the same common frequency and having a constant phase shift. When the amplitudes of both 
oscillators are constant and equal to each other, phase approximation results in a system of Ku- 
ramoto oscillators with distributed-delay coupling, in which case it is possible to find the range of 
admissible collective frequencies and the phase shift analytically. For a uniform delay distribution, 
the system of coupled Kuramoto oscillators exhibits both stable and unstable branches of phase- 
locked solutions, which appear for higher values of the coupling strength as the frequency detuning 
increases. For weak (exponential) delay distribution all the branches of phase-locked solutions are 
unstable, and for the strong (gamma function) delay distribution kernel they are stable for small 
detuning and sufficiently small coupling strength. We have also computed numerically branches 
of phase-locked solutions for uniform and gamma distributions for the full system. For a uniform 
delay distribution, there is co-existence of stable and unstable phase-locked branches, with branches 
being stable for higher values of K as the detuning increases. In the case of gamma distribution, 
there are two branches of phase-locked solutions, which are both unstable, and their frequency 
difference from the average frequency uJ increases with increasing A. 

So far, we have investigated phase-locked solutions in the system of Stuart-Landau oscillators 
with distributed-delay coupling. One possible extension of this work would be the analysis of a 
phase-flip transition (Prasad 2005, Karnatak et al. 2010), where in-phase and anti-phase phase- 
locked branches exchange stability. This phenomenon has been observed in systems with discrete 
time delays, and it has even been observed in the transient dynamics preceding amplitude death. 
Another possibility would be to consider whether coupled systems with delay-distributed coupling 
are able to exhibit other types of phase dynamics and synchronization. One such scenario, which 
is useful in laser applications, is the case when the sum of phases of the two oscillators is constant 
(Erneux &: Glorieux 2010, H.-J. Wiinsche et al. 2005). Phase approximation in this case results in 
a delayed Adler equation, and it would be both theoretically and practically important to consider 
possible solutions of this model for different delay distributions. Other more complex locking 
scenarios have been found in quantum-dot lasers under optical injection (Lingnau et al. 2012, 
Pausch et al. 2012). 
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Appendix 

A convenient way to derive the characteristic equation in the case of gamma distributed delay 
kernel is to use the linear chain trick (MacDonald 1978), which allows one to replace the original 
equation by an equivalent system of {p + 1) ordinary differential equations. To illustrate this, we 
consider a particular case of system (pQ) with a weak delay kernel given by (115H with p = 1, which 
is equivalent to a low-pass filter (Hovel and Scholl 2005): 

g w {u) = ae~ au . (26) 
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(27) 



Introducing new variables 

poo 

Yi{t) = / ae- as Zl {t- s)ds, 
Jo 

roo 

Y 2 (t) = / ae~ as z 2 {t - s)ds, 
Jo 

allows us to rewrite the system ([I]) as follows 

h{t) = (l + i^z^-lzt^fz^ + Ke^^^-zxit)], 

z 2 (t) = (l + iio 2 )z 2 (t)-\z 2 (t)\ 2 z 2 (t) + Ke ie [Y 1 (t)-z 2 (t)}, 
Yi(t) = az-i(t) - aYi{t), 

Y 2 (t) = az 2 {t) -aY 2 (t), 

where the distribution parameter a is related to the mean time delay as a = l/r m . The trivial 
equilibrium z\ = z 2 = of the original system ([1]) corresponds to a steady state z\ = z 2 = Y\ = 
Y 2 = of the modified system (|27p • The characteristic equation for the linearization of system (|2T[) 
near this trivial steady state is given by 

(a + A) 2 (l + iwi - Ke ie - \\ fl + iu 2 - Ke ie - \\ = K 2 a 2 e 2ie , (28) 

which is the same as equation (|17|) , For larger values of p, one would have to introduce a larger 
number of additional variables in a similar fashion (two additional variables for each increase of p 
by l). 
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